#############################################################
############   Additional Balancing ##############
###############      October 10, 2018     ##################
########  Rerun: December 16, 2022 ######

rm(list=ls())
library(coefplot)
library(stargazer)
library(ggplot2)
library(foreign)

#################################################
########### Multiplot Function #################
######   Source: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/



multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
########################################################


#########################################
##### Create  Bandwidths ###############
geo2001=read.csv("~/Dropbox/Personal Research 2017/replications/geocensus2001_nov16.csv")

#Distance to Mysore-Bombay Border
rd10.mb=geo2001[which(geo2001$NEAR_DIST_border1<10000),] #20 km
rd25.mb=geo2001[which(geo2001$NEAR_DIST_border1<25000),] #50 km
rd50.mb=geo2001[which(geo2001$NEAR_DIST_border1<50000),] #100 km
rd100.mb=geo2001[which(geo2001$NEAR_DIST_border1<100000),] #200 km

table(rd10.mb$border1)

#Distance to Hyderabad-Bombay Border
rd10.hb=geo2001[which(geo2001$NEAR_DIST_border2<10000),] #20 km
rd25.hb=geo2001[which(geo2001$NEAR_DIST_border2<25000),] #50 km
rd50.hb=geo2001[which(geo2001$NEAR_DIST_border2<50000),] #100 km
rd100.hb=geo2001[which(geo2001$NEAR_DIST_border2<100000),] #200 km
table(rd10.hb$border2)



#######################################################
#### Plots of space for the balance and unbalance ####
##########     using long and lat    ##################

######borders######
mys=read.dbf("mysborderpts.dbf")
names(mys)
summary(mys)

hyd=read.dbf("hyd_borderpts.dbf")
names(hyd)

#bw 100 km (basically the whole territory)
#mys
g1=ggplot(rd100.mb, aes(x=Longitude, y=Latitude, colour=TOT_POP))+
  geom_point()+
  geom_smooth(data=mys, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g2=ggplot(rd100.mb, aes(x=Longitude, y=Latitude, colour=TOT_SC))+
  geom_point()+
  geom_smooth(data=mys, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g3=ggplot(rd100.mb, aes(x=Longitude, y=Latitude, colour=TOT_ST))+
  geom_point()+
  geom_smooth(data=mys, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g4=ggplot(rd100.mb, aes(x=Longitude, y=Latitude, colour=TerrainRug))+
  geom_point()+
  geom_smooth(data=mys, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

ggplot(rd100.mb, aes(x=Longitude, y=Latitude, colour=Slope))+
  geom_point()+
  geom_smooth(data=mys, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)


multiplot(g1, g2, g3, g4, cols=2)




#hyd
g1=ggplot(rd100.hb, aes(x=Longitude, y=Latitude, colour=TOT_POP))+
  geom_point()+
  geom_smooth(data=hyd, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g2=ggplot(rd100.hb, aes(x=Longitude, y=Latitude, colour=TOT_SC))+
  geom_point()+
  geom_smooth(data=hyd, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g3=ggplot(rd100.hb, aes(x=Longitude, y=Latitude, colour=TOT_ST))+
  geom_point()+
  geom_smooth(data=hyd, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

g4=ggplot(rd100.hb, aes(x=Longitude, y=Latitude, colour=TerrainRug))+
  geom_point()+
  geom_smooth(data=hyd, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)

ggplot(rd100.hb, aes(x=Longitude, y=Latitude, colour=Slope))+
  geom_point()+
  geom_smooth(data=hyd, aes(x=Longitude, y=Latitude), col="green", size=1, se=FALSE)


multiplot(g1, g2, g3, g4, cols=2)



#baseline bw=10km
#mys
g1=ggplot(rd10.mb, aes(x=Longitude, y=Latitude, colour=TOT_POP))+
  geom_point()+
  geom_point(data=mys, aes(x=Longitude, y=Latitude), col="red", size=1)

g2=ggplot(rd10.mb, aes(x=Longitude, y=Latitude, colour=TOT_SC))+
  geom_point()+
  geom_point(data=mys, aes(x=Longitude, y=Latitude), col="red", size=1)

g3=ggplot(rd10.mb, aes(x=Longitude, y=Latitude, colour=TOT_ST))+
  geom_point()+
  geom_point(data=mys, aes(x=Longitude, y=Latitude), col="red", size=1)

g4=ggplot(rd10.mb, aes(x=Longitude, y=Latitude, colour=TerrainRug))+
  geom_point()+
  geom_point(data=mys, aes(x=Longitude, y=Latitude), col="red", size=1)

ggplot(rd10.mb, aes(x=Longitude, y=Latitude, colour=Slope))+
  geom_point()+
  geom_point(data=mys, aes(x=Longitude, y=Latitude), col="red", size=1)


multiplot(g1, g2, g3, g4, cols=2)



#hyd
g1=ggplot(rd10.hb, aes(x=Longitude, y=Latitude, colour=TOT_POP))+
  geom_point()+
  geom_point(data=hyd, aes(x=Longitude, y=Latitude), col="red", size=1)

g2=ggplot(rd10.hb, aes(x=Longitude, y=Latitude, colour=TOT_SC))+
  geom_point()+
  geom_point(data=hyd, aes(x=Longitude, y=Latitude), col="red", size=1)

g3=ggplot(rd10.hb, aes(x=Longitude, y=Latitude, colour=TOT_ST))+
  geom_point()+
  geom_point(data=hyd, aes(x=Longitude, y=Latitude), col="red", size=1)

g4=ggplot(rd10.hb, aes(x=Longitude, y=Latitude, colour=TerrainRug))+
  geom_point()+
  geom_point(data=hyd, aes(x=Longitude, y=Latitude), col="red", size=1)

ggplot(rd10.hb, aes(x=Longitude, y=Latitude, colour=Slope))+
  geom_point()+
  geom_point(data=hyd, aes(x=Longitude, y=Latitude), col="red", size=1)


multiplot(g1, g2, g3, g4, cols=2)





